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We outline a novel approach to develop an in-medium shower Monte-Carlo event generator based 
on the higher-twist formalism of jet modification. By undoing one of the light-cone integrals which 
sets the corresponding light-cone momentum to be equal in the amplitude and the complex con- 
jugate, we introduce an uncertainty in the smaller light-cone momentum component. This allows 
for the generalization of the standard analytic formalism to a Wigner transform like formalism, 
where the non-conjugate large light-cone momentum and position are retained for each parton. 
Jets are generated event-by-event by simulating this Wigner transform kernel. Simple results for 
longitudinal distributions of partons and hadrons from jets propagating through a hot brick of 
strongly interacting matter are presented. Values of the transport coefficient q are dialed to match 
phenomenologically relevant cases. 

PACS numbers: 



I. INTRODUCTION 

With the availability of very high energy jets in 
heavy-ion collisions at the Large Hadron Collider 
(LHC) [IH3], the study of jet modification has moved 
far beyond the suppression of leading hadrons. In 
order to theoretically calculate the large number of 
new measurements such as multi-particle correla- 
tions, jet shapes, flow of energy-momentum in and 
out of jet cones etc., one requires a Monte-Carlo 
based event generator. Even reconstructed single 
particle observables such as the fragmentation func- 
tion of the remnant jets as measured by the CMS 
and ATLAS collaborations [H [5] are, in principle, 
multi-particle effects, where the jet has to be recon- 
structed to some extent, from the parton shower. A 
proper description of all such processes requires a 
Monte-Carlo event generator. 

Monte-Carlo based event generators represent one 
of the most important interfaces between theory 
and experiment in high energy processes, especially 
those involving QCD jets. In vacuum these are very 
well developed [SHE]- So far the development of 
jet modification in the presence of a medium has 
been somewhat varied. This is mostly due to the 
lack of consensus on the right analytic approach to 
jet modification. Prior to the commissioning of the 
LHC, there existed four different analytical energy 
loss approaches [5HP2"] which were equally success- 
ful at describing high transverse momentum (high 
Pt) hadron suppression data from the Relativistic 
Heavy- Ion Collider (RHIC). All of these were based 
on perturbative QCD (pQCD). 

Two of these analytical calculations, enhanced 
with the inclusion of dynamical media, the Armesto- 
Salgado- Wiedemann (ASW) scheme [13] and the 
Arnold-Moore- Yaffe (AMY) scheme [T3] have been 
incorporated within Monte-Carlo event generators. 



The ASW scheme has been incorporated within Q- 
PYTHIA [15], while the AMY scheme has been in- 
corporated within the MARTINI Monte-Carlo event 
generator [16] , Along with these, there are also the 
first principles event generators: JEWEL [I?] [18] 
and YaJEM [TO] [20], which are not completely 
based on any of the analytical models, and in- 
corporate medium effects within vacuum parton 
showers by modifying various matrix elements by 
hand. This was also the methodology followed 
by the very first jet modification event genera- 
tor PYQUEN [21 , where the energy loss kernel 
was loosely based on the Baier-Dokshitzer-Mueller- 
Peigne-Schiff (BDMPS) approach [22]. 

The generic methodology in all these cases has 
been to take the standard vacuum event generator 
and either add a component that represents medium 
effects, after the full vacuum shower has been gen- 
erated, or to enhance the vacuum shower algorithm 
itself such that the resulting shower is modified, i.e., 
vacuum and medium induced radiation is accounted 
for simultaneously. Besides the variety of techni- 
cal issues that must be surmounted in the construc- 
tion of an in-medium jet Monte-Carlo event gener- 
ator, two fundamental problems need to be dealt 
with, which have no analog in vacuum: The intro- 
duction of space-time structure (and fluctuations in 
space-time structure) in the shower, and a mecha- 
nism of hadronization where partons from the hot 
dense medium will coalesce with partons from the 
jet to form hadrons. None of the of the above men- 
tioned schemes addresses either of these two effects 
in a satisfactory manner. 

In this paper, we attempt to solve the first of these 
problems: the incorporation of a theoretically con- 
sistent space-time structure, as well as fluctuations 
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in space-time structure, 1 within the event generator 
framework. The Monte-Carlo event generator will 
be based on the Higher- Twist scheme, and this pa- 
per will contain the first set of results within such a 
framework. The event generator will be referred to 
as the Modular- All-Twist- Transverse-and-Elastic- 
scattering-induced-Radiation (MATTER++) event 
generator. The '++' simply indicates the choice of 
programming language. The remaining terms in the 
moniker refer to the underlying theory which may 
be found in Refs. [H|25]. 

The remaining paper will be organized as follows: 
Sec. II will contain a terse review of final state event 
generation for jets in vacuum. Sec. Ill will discuss 
how the light-cone location of a split, which is nor- 
mally integrated out, may be reintroduced into the 
calculation of jet modification. In Sec. IV, we will 
outline the basic components of our Higher- Twist 
based event generator. In Sec. V, we will present 
results of our numerical calculations of jet modifica- 
tion in a static medium of finite length. Concluding 
discussions are presented in Sec. VI. 



II. FINAL STATE VACUUM EVENT 
GENERATION 

In this section we will review the simulation of 
parton showers in vacuum. In order to illustrate 
the reintroduction of space-time in jet shower de- 
velopment, we will set up the process within the 
framework of Deep-Inelastic Scattering (DIS). We 
will take a factorized approach, valid at high ener- 
gies, in the presence of a hard scale. We will ignore 
the development of the initial state space- like shower 
and only consider the final state time-like shower. 

A hard lepton with incoming momentum Li scat- 
ters off a proton with incoming momentum P and 
exits with an outgoing momentum L. The proton 
is shattered into a large number of final hadronic 
states. The inclusive differential cross-section for 
an arbitrary process with hadrons with momenta 
Pi . . . pn , in the single photon approximation, may 
be expressed as, 
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where, is the leptonic tensor, all factors which 
represent the interaction of the incoming and outgo- 
ing lepton with the exchanged photon are contained 
within . The hadronic tensor W^ v contains the 



1 Note: fluctuations in space-time structure have been phe- 
nomenologically introduced in YaJEM 23 . 



entire set of hadronic interactions including the in- 
teraction of the proton with the photon. The Man- 
delstam variable s = (xbP + q) 2 ■ The Bjorken vari- 
able xb = Q 2 /(2P ■ q). We choose a frame where 
the momentum of the photon is given in light-cone 
coordinates as, 
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often referred to as the Breit frame. 

The simplest instance of a hadronic final state is 
a single quark in the final state. In this case the 
hadronic tensor has a simple factorized form, 
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In the equation above, Q q represents the charge of 
the quark in units of the electron's charge, f q rep- 
resents the light-cone distributions of quarks inside 
the proton. 

The equation above is the integral hadronic ten- 
sor. The differential hadronic tensor as a function 
of the three dimensional momentum of the outgoing 
quark is given as, 



d 3 W^ 
d 3 l 



= WrS 2 (lx)S(r ~q-), 
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where, the outgoing quark carries the (— )-lightcone 
momentum of the incident photon. The equation 
above represents a simple expression, where the pro- 
duction cross section of the hard parton, contained 
within Wq V ', is factorized from the final state distri- 
bution, which in this case is simply the three dimen- 
sional delta function S 3 (l — q). In what follows, we 
will not discuss the production process of the hard 
parton and only focus on the final state propaga- 
tion. In the case of a parton produced close to its 
mass shell, there is scarce little beyond the 3-D delta 
function, that may be calculated with pQCD. 

For the case of hard parton, produced off its mass 
shell (time like for a final state parton), one has to 
include multiple final state emissions. In this process 
of emission, the virtuality of the hard parton is grad- 
ually reduced to a low enough scale where hadroniza- 
tion can begin to set in. While the hadronization 
process requires phenomenological modeling, emis- 
sions which lead to the reduction of the virtuality 
from Q 2 to a minimum hard scale where pQCD is 
applicable can be described by constructing and sim- 
ulating a Sudakov form factor. For the case of a 
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quark in vacuum, the Sudakov form factor is given 
as, 



S(Ql,Q 2 )=cxp 



Q 2 l-Qa/Q 

dyPqg(y) 
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In the equation above, P qg {y) is the usual LO split- 
ting function for a quark to radiate a gluon which 
carries a forward momentum of (1 — y)q~ , and have 
its light-cone momentum degraded to yq~ . The 
quantity /i 2 represents the highest allowed virtuality 
of the parton undergoing the split and is also the 
running scale in the process. The Sudakov form fac- 
tor represents the probability that there is no emis- 
sion, between the scales Qq to Q 2 . 

The procedure to generate a shower in vacuum is 
now completely straightforward, and has been de- 
scribed in several texts [27] • One samples the Su- 
dakov form factor and obtains the highest virtuality 
that the remaining process may possess. This is fol- 
lowed by sampling the splitting function, to obtain 
a fraction y which represents the fraction of light- 
cone momentum left in one of the partons after the 
split. This process is then iterated for each virtual 
parton to generate the full shower. The cascading 
process in any part of the shower is stopped when 
the virtuality of that part reaches Q 2 ,. 

In the process of shower development in the vac- 
uum, the entire cascade process may be carried out 
in momentum space. Note that space-time does not 
enter in the formulation of the Sudakov form factor 
in vacuum [Eq. In what follows, we will con- 

sider the modification of this process in the presence 
of a medium; here the location of a hard parton will 
need to be retained at all stages of the calculation. 
The modification of the shower Monte-Carlo algo- 
rithm will be discussed in Sect. IV. 



III. THE RE-INTRODUCTION OF 
SPACE-TIME IN MONTE-CARLO EVENT 
GENERATION 



the gluon field generated in a medium. In this limit, 
the multiple scattering expression, for a quark prop- 
agating through the medium without emission, can 
be expressed using almost classical arguments as the 
gradual linear broadening with distance travelled. 
At some point in its path, the quark will radiate a 
gluon and these two partons will, beyond the for- 
mation length, propagate as independent particles 
multiply scattering off the gluon field. At the split, 
the two partons will separate from each other with 
a large and opposing transverse momenta l±. The 
accrued transverse momentum k± adds incoherently 
to this large momentum 

To carry out the program of broadening described 
above, the crucial piece of information required is 
the (light-cone) location, beyond which the two par- 
tons that have split from the parent parton, may be 
considered to scatter independently. The reader will 
note that while the location of a split, on average, 
will be different between jets in vacuum versus those 
in medium, this remains a well defined question even 
for the case of jets in vacuum. We will first attempt 
to answer the question of the location of splits for 
jets in vacuum. 
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FIG. 1: The diagram for the hadronic tensor of DIS on 
a proton with an outgoing quark and a radiated gluon. 



In the preceding section, the pQCD based vac- 
uum shower Monte-Carlo formalism was reviewed. 
In this well tested set up, no information about the 
space-time location of the shower is retained. In this 
section we will make the most minimal modification 
to the formalism to include this information. In the 
subsequent section, we will describe the modification 
to the propagation of non-emitting partons, and that 
to the Sudakov form factor due to scattering in the 
medium. 

Consider, for simplicity, the process of high en- 
ergy DIS with a near on-shell quark, scattering off 



Consider the process of semi-inclusive DIS on a 
proton, with a hard quark and a radiated gluon in 
the final state (illustrated in Fig. [T]). The hadronic 
tensor may be expressed as, 

Wr = J d 4 yoWM7^7>(0)|P) 

= J dV,TM^7^7"]^(lfo)0(»b). (6) 

The location of the quark in the complex conjugate 
is yo and its location in the amplitude is set as the 
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origin. The factor F(y ) represent the expectation 
of operators that constitute the initial state, while 
O(yo) constitutes the expectation of operators in the 
final state. 



In DIS, the momentum components of the hard 
photon q~ and q + = Q 2 /(2q~) are experimen- 
tal measurements. In most analytical calculations, 
any uncertainty in these quantities is averaged over. 
Since, no attempt to determine the space-time loca- 
tion of the hard scattering is made, the uncertainty 
in these quantities is minimal. However, this is not 
the case for p£ = xbP + ■ The value of Xb is an ap- 
proximate average value. Hard quark jets are pro- 
duced far off their mass shell. The cause of this 
off-shellness is the uncertainty in p^ . The formation 
length of any emission in the jet's shower develop- 
ment depends on the virtuality of the emitting par- 
ton; more precisely, it depends on the uncertainty 
in the virtuality, limited by the virtuality itself. To 
elucidate the space-time structure of the jet created 
by this hard interaction, the uncertainty in the var- 
ious momentum components has to be taken into 
account. 



No doubt, all other energy- momentum and space- 
time factors will exhibit uncertainties. However, as 
Pq represents the largest scale in the problem, the 
uncertainty in this quantity has the maximal effect 
on the space time structure of the developing jet. Let 
us denote the momentum of the final state quark as 
p = [p + ,p~ ,p±\- Quark jets produced in DIS have 
a large (— )-light-cone momentum component and a 
small (+)-component. As a result, the space-time 
profile of the jet is concentrated near the x + = 
light-cone. In this first attempt, we will not probe 
the deviation of the jet away from the strict light- 
cone. As a result, we will ignore the uncertainty in 
p~ (or pj_) and only focus on the uncertainly in p + . 
Note that the uncertainty in p^ directly translates 
into an uncertainty in p + . In a future effort both 
uncertainties will have to be considered. An uncer- 
tainty in p is most straightforwardly introduced by 
setting the p vector to be different in the amplitude 
and complex conjugate. In the end, the offset in all 
components of p except in p + will be ignored. 

In the presence of this uncertainty in p, the initial 
state operator expectation F(y ) is unchanged: 



The final state operator expectation becomes 



f(vo) 



{pm )\m\p), 



(7) 



O 



Tr 



2_ 

2 



O 



(8) 



d 4 zdV d% ^ dVo 



x Tr 



(2tt) 4 (2tt) 4 (2tt) 4 (2tt) 4 

7" ~«(^o+ v ' 



2 {p + q) 2 -ie 



x G a p(l)2Tr5(l 2 )(-i^)-^ 



(p'o + q) 2 + it 

x gil-Vo e -»(Po+«)-(vo— 2) e -U-{z—z') e -il r (z-z') 

x e- l( - p '° +qyz ' g 2 . 

As pointed out in the diagram in Fig. [TJ the incom- 
ing quark momentum in the amplitude and complex 
conjugate (p' and po) are mismatched. Similarly, 
the location of the emission of the gluon is denoted 
as z' in the amplitude and z in the complex conju- 
gate. Note that they do not have to be at the same 
point. 

The phase factors and the z, z' integration may be 
simplified as follows, 
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In the equation above, po = (po +J 5 o)/2 and Apo = 
Po ~ Po- Similarly the two positions are defined as 
z = (z+z')/2 and Az = z—z'. Carrying out the Az~ 
integral we obtain the usual relation that p + +q + = 
1+ + 1+. 

In the limit that Pq and p' Q are vanishingly small, 
the matrix element may be simplified as, 
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In the equation above, p°y° is expressed as Ap°yo + 
PoAy . Since no part of Eq. ^ depends on p^ 
or Pq , these have been integrated out to obtain 
5{z + )8{z' + ), localizing the process on the light cone. 
Note, this would be modified if Aq~ had been re- 
tained. This leaves the transverse momentum inte- 
grals. Carrying out the z± integrals, sets po_L = Po±- 
One may also integrate out the Apo_L which yields 
a S- function, that may be used to set z± = j/o_i_- 

The final integration that is left is that over z~. 
In order to carry it out, one needs to evaluate the 
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ApQ distribution. To compute this exactly, one has 
to carry out the integral, 
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Note that in Eq. ([6]), there is no integration over the 
mean distance yo- This is because, we have set the 
location of the scattering in the amplitude to be the 
origin. But in reality, all locations in the nucleon 
are not equivalent and one should average over this 
variable. The result of such an averaging may be ex- 
pressed as the usual parton distribution function for 
Pq times a distribution for the offset ApJ. For this 
first attempt, we will assume a Gaussian uncertainty 
for the A-Pq distribution: 
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p(Ap + ) = 



V2~7 



(12) 



The width depends on the scale of the process (Q 2 ): 
Processes with higher Q 2 or rather higher q + have 
larger widths. We will set the width using a boot- 
strap method where, given a width, we first calcu- 
late the distribution of locations z~ where the split 
occurs, and then average over these locations to ob- 
tain the mean distance of the split. The distribution 
of z~ is obtained by Fourier transforming the dis- 
tribution for Ap + , which will also yield a Gaussian 
distribution for z~ (Note: the distribution of Ap + 
is identical to the distribution of ApQ ). The width 
is then dialed so that this mean distance equals the 
formation time for a given q~ and Q 2 . 

As a result of these simplifications, the produced 
off-shell quark has a momentum of p = [p~~ ,p + , 0, 0]. 
We obtain the forward light-cone momentum p~ = 
q~ and a virtuality fi 2 — 2p + p~ = 2p$q~, where 
Pq is the momentum of the incoming quark. The 
gluon radiated from this quark will have a light-cone 
formation time (equal to the light-cone formation 
length for a hard parton), 



pf 



(13) 



Thus, the Ap distribution will have a width of 
a = 2p + /tt = 2pQ /n. This will ensure that on in- 
tegrating the Gaussian distribution from a distance 
z~ = to oo, the mean light-cone distance will be 
equal to the light-cone formation time in Eq. (13). 
To obtain a distribution of formation times, we sim- 
ply have to sample this Gaussian distribution for 
positive values of the light-cone location of the split. 

While motivated for the case of the first emis- 
sion in DIS, the set up may be used for the case of 



multiple emission: One simply determines the large 
light-cone momentum p~ and the off-shellness /x 2 , 
or rather the smaller light-cone momentum p + = 
/i 2 /(2p~), and then computes the formation time. 
Assuming that the uncertainly in p + is of the order 
of (but smaller than) p + itself, we obtain a distri- 
bution in formation times by Fourier transforming 



Eq. (12) 
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with a = \/2p + /y/TT. The inclusion of this distri- 
bution in x~ , in combination with the distribution 
in p + which will be obtained by sampling the Su- 
dakov form factor convert the shower formalism into 
a Wigner Transform like formalism where the two 
non-conjugate components of space-time and mo- 
mentum are retained. These will be used to deter- 
mine the location of the split of a particular parton. 

In the set up described above, to obtain a particu- 
lar splitting length, one needs to know the two com- 
ponents p~ and p + . In the subsequent section we 
will describe how to obtain these momentum com- 
ponents p~ and p + from the calculation of the Su- 
dakov form factor. The Sudakov form factor will be 
determined for the case in vacuum and in the case 
of a medium. The effect of the medium will enter 
only in the determination of the distribution of p + 
and p~ . Once p + and p~ have been determined, 
the distribution in Eq. ( 14 ) can be used to find the 



splitting distance for a particular event. 



IV. FINAL STATE IN-MEDIUM EVENT 
GENERATION 

In the preceding section, the means to calculate 
the distribution of splitting lengths was outlined. 
We have refrained from denoting the lengths ob- 
tained from sampling Eq. ( 14 ) as formation lengths 
(or formation times). The mean of the splitting 
lengths will be denoted as the formation length (or 
formation time) . In the calculation of the distribu- 
tion of splitting lengths, the input is composed of the 
light-cone momentum components of the hard par- 
ton which will undergo the split. In this section we 
will outline how these components are determined. 

We return to the case of a hard quark produced 
in DIS on a nucleon, and simply replace the nucleon 
with a large nucleus. The photon hits a quark in one 
of the nucleons and then sends it back through the 
nucleus where its evolution to a jet is modified by 
the presence of a nuclear medium. One may divide 
the modification into two parts, one which modifies 
the transverse momentum distribution of a propa- 
gating parton without inducing it to radiate and an- 
other which changes the radiation probability of the 
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virtual parton. As in the case of the analytic calcu- 
lation [Ml [25] , we will compute these separately and 
these within one calculation. 



lation [24J |28] , we will compute these separately 
then incorporate these within one calculation. 

The modification of the transverse momentum dis- 
tribution, of the outgoing quark, without emission, 
in a medium, is trivial: As pointed out in [25 , 28 , 29 , 
this may be easily expressed by simply replacing the 
(5-function by position dependent Gaussians. In this 
first attempt to convert the Higher- Twist approach 
to a Monte-Carlo, we will ignore elastic drag. 

In the following, the production of the hard par- 
ton will no longer be considered and the focus will 
lie solely on the fate of the final state parton. Trans- 
verse momentum broadening is achieved by sampling 
a unitary Gaussian distribution. The parton enters 
the medium with a narrow distribution in trans- 
verse momentum [<5(Zj_)], which is then broadened 
to a two dimensional unitary Gaussian distribution, 
whose width depends on the length traversed by the 
parton, 
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In the equation above, D is the transverse momen- 
tum diffusion coefficient, the well known transport 
coefficient q = AD. The light cone length L~ is the 
light-cone length traversed in their medium. In the 
case that the medium is not uniform, one may sim- 
ply replace DL~ — > J dL~ D(L~). 

To determine the lengths traversed in the medium, 
one requires the locations of the the various splits 
that a jet undergoes as it propagates through the 
medium. To be specific, let us consider a single 
hard quark propagating through a dense medium un- 
dergoing multiple hard scattering, and both vacuum 
and medium induced radiation. As mentioned in the 
preceding section, to determine the location of the 
splits, one needs the two light-cone components of 
the virtual quark. The (— )-component p~ will be 
considered to be large and the (-l-)-component p + 
will be considered to be small (but still much larger 
than Aqcd)- The virtuality Q 2 = 2p + p~ . 

In order to obtain the distribution of leading 
hadrons emanating from the fragmentation of this 
hard quark, one would compute the medium mod- 
ified fragmentation function [30] . This requires the 
solution of the Dokshitzer-Gribov-Lipatov-Altarelli- 
Parisi (DGLAP) evolution equation I3"TH3"3] in 
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In the equation above D^(z, Q 2 ,p~)\^ is the medium 
modified fragmentation function for a hadron to to 
fragment with light-cone momentum zp~ (0 < z < 
1), from a hard quark with a light-cone momentum 
p~ , and virtuality Q 2 . For this fragmentation func- 
tion, the quark commences propagation at Q and 
exits at £/■ The factor P(y) represents the un- 
renormalized splitting function, and y represents the 
fraction of momentum which remains in the quark 
after a gluon has been radiated. The leading twist 
contribution to the multiple scattering, single emis- 
sion kernel K, is given as 
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The full evolution equation for the medium mod- 
ified fragmentation functions will include contribu- 
tions from gluon fragmentation functions which have 
a similar in-medium evolution. 

We should point out at this stage, that in the 
Higher- Twist approach, multiple emissions are or- 
dered in transverse momentum as they are in vac- 
uum DGLAP emission. Multiple scattering is sup- 
pressed by powers of the hard scale and enhanced 
by the length of the medium traversed by a par- 
ton. Such ordered multiple emission in the pres- 
ence of multiple scattering is only considered when 
qL <C Q 2 (L is the length traversed by parton be- 
fore splitting) i.e., when the multiple soft scattering 
does not produce a considerable affect on the trans- 
verse momentum distribution between the partons 
post split (or the virtuality distribution of the par- 
ent parton). Once Q 2 becomes too small, Higher- 
Twist based calculations have to be abandoned. We 
should point out that beyond this limit, the ordering 
of multiple emissions is, as yet, unsettled |34| . 

Since the medium modified evolution equations 
are valid in a similar range of momenta as the vac- 
uum emission process, we may set up a Sudakov form 
factor for the Higher- Twist based multiple scattering 
induced emission process. In this case, the medium 
modified Sudakov form factor for a quark with light- 
cone momentum p~ , with a maximum virtuality Q 2 , 
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propagating through a dense medium, starting at 
the light-cone location C~ is given as, 



5 cr (Qg,Q 2 ) = exp 
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In the equation above, the vacuum portion of the ar- 
gument of the exponent is identical to Eq. ([5]). The 
medium portion is integrated from the location C~ 
to the location £r +t~ . The choice of the distance r 
is significant. There are two ways to deal with this: 
One method is the use the sampled value of Q 2 and 
the known value of forward light-cone momentum 
p~ (which is equal to q~ in this case) and substi- 
tute p + = Q 2 / (2p~) in Eq. (14). One can then sam- 



ple this distribution and obtain a particular splitting 
length x~ . This length may then be substituted in 
Eq. ( 18 1 to compute the Sudakov factor. The down- 



side of this method is that one will always get a ran- 
dom length x~ for a given step in the Monte-Carlo 
simulation [averaging over many such samplings of 
Eq. (14 1 will yield the mean formation lenght]. As a 



result the Sudakov from factor will not be a mono- 
tonic function of Q 2 which will obviously slow down 
the computational speed of any such calculation. 

An alternative is to use a mean length, which in 
this case would be the formation length tJ based on 
p~ and the chosen virtuality Q 2 . One may be criti- 
cal about this mean choice, especially since one is us- 
ing the mean of the vacuum distribution of splitting 
lengths as the length of the medium modified part 
of the calculation as well. However, this choice of 
r~ = tJ may also be justified based on the medium 
modified calculation. In Fig. [2j we plot the medium 
modified single emission kernel. The red solid line is 
a plot of K from Eq. ( 17), modulo the constant fac- 
tors 2q/Q 2 , as a function of light-cone length x~ (in 
units of the light-cone formation length t7). This 
is referred to as the GW line as it was first derived 
for the case of single scattering induced radiation by 
Guo and Wang 0[3S]. 

The GW formula has been derived for the case of 
short lengths in the medium and has the property 
of continuously growing with distance. One obtains 
this formula by ignoring a series of phase factors 
which are small for short distances [3B]. However, 
the GW formula yields no obvious point to stop the 
length integration. This is resolved by using the for- 
mula derived by Aurenche-Zakharov-Zaraket [571155] 
(referred to as AZZ and plotted as the dotted line in 
Fig. [2]). This formula contains all the phase factors 
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FIG. 2: A comparison between the behavior of the 
Guo-Wang (GW) and the Aurenche-Zakharov-Zaraket 
(AZZ) radiation kernels as a function of length traversed 
by parton in a medium. 



that were ignored by Guo and Wang. It matches 
with the results of GW for distances up to the for- 
mation length. Beyond the formation length, it de- 
viates from the GW results, turns over and starts 
dropping down to zero. Beyond twice the forma- 
tion length, the AZZ result becomes negative and 
then oscillates continuously up to infinity. The in- 
terpretation of this result, is that the single emission 
formula (or perturbation theory up to this order) is 
only valid up to a distance of the order of the for- 
mation length. Beyond this length, one must use 
multiple emissions. Based on these considerations, 
we chose r~ = tJ and the GW formula for the ker- 



nel in Eq. (18). One may also chose r~ = 2r7 and 
the AZZ formula to retain the full range of posi- 
tive AZZ contributions. For a static homogeneous 
medium this can be easily equated with the previ- 
ous prescription by rescaling q. 

With the Sudakov factor, now completely defined, 
it may be sampled to obtain the virtuality of a par- 
ticular parton. The large light-cone momentum p~ 
and the virtuality Q 2 can now be used to sample 
the splitting length distribution in Eq. (14) to ob- 



tain an event-by-event distance where the parton 
will undergo a split. The remaining steps are identi- 
cal to those in vacuum: one computes the splitting 
fraction y by sampling the splitting function. This 
gives two partons with light-cone momentum yp~ 
and (1 — y)p~ , and transverse momentum l± and 
— Ixj where l± — Q 2 y(l — y). The maximal virtu- 
ality of the two produced partons is Q ; the actual 
virtuality will be obtained by again sampling the Su- 
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dakov factor for each of these separately. 

With the setup described above, we can now gen- 
erate a parton shower. With the knowledge of the 
location of the splits, the distance traveled by a par- 
ton between splits is known. This allows for a calcu- 
lation of the transverse broadening by scattering to 
be carried out using Eq. ( 15 1. In the subsequent sec- 



tions, we will present simple numerical calculations 
of jet modification in a dense extended medium. 



V. NUMERICAL SIMULATIONS 

In the preceding sections, the formalism to carry 
out event generation of a parton shower propagat- 
ing through the medium in MATTERH — h was set up. 
In this section, simple numerical simulations will be 
presented. In order to study the various facets of jet 
propagation in a dense medium, independent from 
the evolution of the medium itself, in this first effort, 
we will only consider jet propagation in a finite sized 
static medium held at a fixed temperature. This set 
up is often referred to as the brick [39 . Event gener- 
ation has been studied within a brick format by sev- 
eral other authors [TTJ HO] • In the remainder of this 
section, we outline 3 separate calculations: The cal- 
culation of the distribution of surviving quarks and 
gluons emanating from a single hard quark or gluon, 
the medium modified fragmentation function, and a 
medium modified fragmentation function based on 
the surviving energy of a hard jet (motivated by re- 
cent CMS and ATLAS results). 

The first test of our new Monte-Carlo (MC) for- 
malism will be a comparison of the distribution of 
partons (with a minimal virtuality of Qo — 1 GeV) , 
as a function of the longitudinal light-cone fraction 
z = Ph/p~, in the vacuum with standard analytic 
calculations. In Fig. |3j we consider the case of a 
hard quark with an initial q~ = 100 GeV and a 
maximal virtuality of Q = 50 GeV. This leads to 
a shower of partons in the vacuum. The virtuality 
drops with each step of the shower until a parton 
reaches Q = 1 GeV. The parton stops branching 
at this point. The final distribution of quarks and 
gluons represents this collection of partons with a 
Q = 1 GeV. The hollow circles represent the distri- 
bution of quarks, while the hollow squares represents 
the distribution of gluons. Results are presented for 
100,000 events. 

Comparison with analytic DGLAP calculations is 
not straightforward. One starts with a quark-to- 
quark and gluon-to-quark fragmentation function at 
the lower scale of Q 2 = 1 GeV 2 , 



D q ^. q (z, 1 GeV ) = 8(1 -z) 
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FIG. 3: (Color online) A comparison between the vac- 
uum MC calculation and a DGLAP calculation of the 
distribution of partons with Q = 1 GeV, as a function 
of the light-cone momentum fraction z. The originating 
parton is a 100 GeV quark. MC results are for 100,000 
events. See text for details. 



This is then evolved up to Q = 50 GeV to obtain the 
LVhj(z,2500 GeV 2 ) and D g ^ q (z, 2500 GeV 2 ). The 
quark-to-quark fragmentation function is plotted as 
the solid blue line in Fig. [3| One then performs 
another calculation where the starting point at Q = 
1 GeV is given as, 



D 



q^ g (z, 1 GeV ) = 0, 



D g ^ g (z,l GeV ) 



8(1 -z). 



(20) 



These are then evolved up to Q 2 = 2500 GeV 2 and 
the quark-to-gluon fragmentation function at this 
upper scale is plotted as the dashed magenta line in 
Fig. [3] In Fig. |4j we reverse the order in the input 
fragmentation functions defined above to consider 
the case of a gluon jet with a Q = 50 GeV, shower- 
ing in vacuum, and terminating with a distribution 
of quarks and gluons at Q = 1 GeV. 

As the reader would have guessed, the analytic 
calculations are hindered by the fact that the 8- 
function is not an acceptable input to the DGLAP 
evolution equations, as it is singular at z = 1. This 
singularity has to be regulated by choosing a particu- 
lar representation of the <5-function. As was noted in 
Ref. [39], only one of the standard regulators where 



D g ^ q (z, 1 GeV 2 



(19) 



5(x) = lim - V \x\ < e, 
8(x) = limO V |x| > e, 



(21) 
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FIG. 4: (Color online) Same as Fig. [3] except that the 
originating parton is a 100 GeV gluon. MC results are 
for 100,000 events. See text for details. 



has yielded convergent results as e — ¥ 0. In spite of 
such issues, the agreement between the MC simula- 
tions and the DGLAP calculations is rather good. 
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FIG. 5: (Color online) A comparison between the 
in-medium MC calculation and an in-medium DGLAP 
calculation. The medium has a length of 5 fm and a 
q — 1 GeV 2 /fm. The originating parton is a 100 GeV 
quark. MC results are for 100,000 events. 

We next present a comparison between the MC 
calculations for a hard jet propagating through a 



dense medium and an in-medium DGLAP calcula- 
tion. These are presented for a parton with an ini- 
tial energy of 100 GeV, a Q = 50 GeV, propagating 
through a medium of length 5 fm with a homoge- 
neous q = 1 GeV 2 /fm. Fig. [5] represents the case 
of the initial parton being a quark, and Fig. [6] rep- 
resents the case of the initial parton being a gluon. 
The comparison is not as good. This is mostly due 
to the approximations introduced in the in-medium 
DGLAP calculations, primary among which is the 
neglect of space-time dependence of the fragmenta- 
tion function [30] . The result of this neglect is that 
the DGLAP evolutions equations show a higher re- 
sult for z — > 1 and are smaller than the MC calcula- 
tions over most of the range of z. 
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FIG. 6: (Color online) A comparison between the 
in-medium MC calculation and an in-medium DGLAP 
calculation. The medium has a length of 5 fm and a 
q = 1 GeV 2 /fm. The originating parton is a 100 GeV 
gluon. MC results are for 100,000 events. 

As we mentioned above, the DGLAP evolution 
equations show regulator dependence when the in- 
put is a (^-function. To remove this defect in the 
comparison, we fold both the MC calculations and 
the DGLAP evolutions equations with an exper- 
imentally measured hadronic fragmentation func- 
tion. In the case of the DGLAP evolution equations, 
this is straightforwardly done by replacing the 5- 
function with a fragmentation function (in this par- 
ticular case, we use a ir° fragmentation function) at 
the lower scale of Qq — 1 GeV 2 . In the interest 
of simplicity, we use the Binnewies-Kniehl-Kramer 
(BKK) fragmentation function [5T] , which has a sim- 
ple analytic expression. This is then evolved up 
to the scale of Q 2 = 3 GeV 2 . For the case of the 
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Monte-Carlo calculation, there is no straightforward 
means to fold a fragmentation function: one needs 
an hadronization model such as string breaking in 
PYTHIA [6l|42] or cluster fragmentation as in HER- 
WIG [7JI33]. Here we use a very simple ansatz: We 
allow the partons to first shower in the medium and 
then continue to shower in the vacuum until they 
reach a Qq = 1 GeV 2 . We then convolute the BKK 
fragmentation function at a Qq = 1 GeV 2 with the 
longitudinal distribution of final partons (as plotted 
in Figs 3][6). Partons that reach the lower scale of 
Qo = 1 GeV 2 while still deep in the medium (more 
than 1 fm from the exit surface) are excised from 
the shower. This hadronization set up should be 
accurate for the z — > 1 (hard) partons which are ex- 
pected to hadronize independently and become less 
accurate for soft partons, where one expects several 
partons to participate collectively in the hadroniza- 
tion process. The results of this comparison are 
plotted in Fig. [7j where we present the ratio of the 
medium modified fragmentation function to the vac- 
uum fragmentation function of 7r°'s, from a hard 
initial quark. In the medium modified calculation, 
we have chosen a medium of length 4 fm with a 
q = 0.2 GeV 2 /fm. This is comparable to the mea- 
sured ratio of the yields of leading hadrons observed 
in the case of DIS on a large nucleus to that in a pro- 
ton, as seen by the HERMES experiment |44 ]. The 
parameters used are somewhat close to the case of 
the Krypton nucleus [3U] . Once again, the agreement 
between the two calculations is quite good. Monte- 
Carlo calculations are carried out for 100000 events. 
There are error bars associated with both the vac- 
uum and medium modified MC simulation, which 
have not been presented. 

Having tested the Monte-Carlo simulation against 
analytical DGLAP calculations, we now proceed to 
calculate quantities that may only be accessed via 
a Monte-Carlo event generator. In Fig. [7J the frag- 
mentation functions are defined in the traditional 
way, where the light-cone fraction z = jp~ is de- 
fined with respect to the original energy of the hard 
parton p~ before it enters the medium. Recently, 
the CMS and ATLAS collaborations have defined 
the fragmentation functions in another way [H [5]. 
The jet with degraded energy is reconstructed and 
one then calculates the fragmentation function with 
the degraded total energy p'~ , i.e., z' = p^/p'~- By 
degraded jet energy, we mean that some of the en- 
ergy of the jet has been lost in the medium and is not 
regained during reconstruction. In all these calcula- 
tions, partons whose virtuality has dropped below 
Q 2 = 1 GeV 2 while still deep inside the medium 
(greater that a fermi away from the exit surface) are 
excised from the shower. This is the case in both 
Fig. [7] and also in the calculations that will be de- 
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FIG. 7: (Color online) A comparison of the ratio of 
medium modified to vacuum fragmentation functions, as 
calculated by the Monte-Carlo event generator and using 
DGLAP evolutions equations. Both calculations take, 
as input, vacuum BKK fragmentation functions at the 
lower scale of Qq — 1 GeV 2 . 



scribed in Fig. [8] 

One notes that the fragmentation functions ob- 
tained by the CMS and ATLAS collaborations are 
quite different (even qualitatively) from those mea- 
sured by the HERMES collaboration, in DIS on a 
large nucleus. For the case of HERMES, the frag- 
mentation functions are defined with respect to the 
original energy of the jet before it enters the medium 
and have a behavior similar to that in Fig.[7j i.e., the 
ratio of the medium modified fragmentation func- 
tion to the vacuum fragmentation function is mono- 
tonically decreasing with increasing z. However, 
the fragmentation functions measured by CMS and 
ATLAS show an initial drop and then a rise, with 
increasing z, in the ratio of the medium modified 
fragmentation function to the vacuum fragmenta- 
tion function. We will demonstrate that this sort 
of non-monotonic behavior is typical of cases where 
the entire energy of the hard jet is not reconstructed, 
i.e., cases where some of the energy is lost in the 
medium. 

In Fig. [HJ we plot various versions of the ratio 
of the in-medium fragmentation functions to vac- 
uum fragmentation functions. In all cases, for the 
medium modified fragmentation function, we have 
an originating quark jet or a gluon jet with an orig- 
inal E = 100 GeV, and a maximum Q = 100 GeV. 
The maroon squares and green crosses represent the 
case where the fragmentation functions are defined 
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FIG. 8: (Color online) Calculations of the ratio of the 
medium modified fragmentation function to the vacuum 
fragmentation functions for a quark jet and a gluon jet 
with an initial energy of 100 GeV, propagating through 
a medium of length 5 fm, with a q = 1 GeV 2 /fm. The 
green crosses and maroon squares represent fragmenta- 
tion functions constructed using the full initial energy 
of the quark and gluon jet respectively. The remaining 
plots represent fragmentation functions constructed us- 
ing the final degraded energy of the jet, which results 
from the non-detection of partons with energies below 
2 GeV and 4 GeV. See text for datails. 



with respect to the original jet energy, as in the case 
of HERMES. These represent the standard defini- 
tion of the fragmentation function (with z! = z). 
and the ratio of medium modified to vacuum shows 
a monotonic decrease with increase in z. The red di- 
amonds and green up-triangles represent quark and 
gluon jets where all partons in the shower with a 
forward light cone momentum below _p~ in = 2 GeV 
(i.e. E < y/2 GeV) are excluded. The jet is recon- 
structed from the ensemble of surviving partons and 
has a lower light-cone momentum p'~ . The light 
cone fraction z' — p~^jp'~ is defined with respect 
to this lower momentum. This automatically moves 
particles with a p~^ = zp~ to a higher value of the 
variable z. Since the fragmentation functions are 
steeply falling with z, this leads to an enhancement 
at larger z. The enhancement is larger for the case of 
gluons as the gluon fragmentation function is more 
steeply falling than the quark fragmentation func- 
tion. Increasing the threshold of partons which sur- 
vive to p~ nin = 4 GeV, leads to a greater shift of 
final hadrons to higher z and results in more en- 
hancement. Similar conclusions were also reached 



by a more sophisticated analysis within the YaJEM 
event generator in Ref. 20J . 

The reader will note that while qualitatively sim- 
ilar, these ratios of fragmentation functions do not 
show quantitatively the same behavior as that ob- 
served by CMS and ATLAS. We would point out 
that these calculations are performed for a static 
brick of finite length, all medium dynamics have 
been ignored. Also one does not observe quark and 
gluon fragmentation functions separately, but rather 
a mixture of the two where the amount of admixture 
depends on the momentum of the detected hadrons. 
In future publications, these Monte-Carlo showers 
will be propagated through a dynamical medium to 
calculate the modification of hard jets. 



VI. CONCLUSIONS 

In this paper, we have presented the first Higher- 
Twist based Monte-Carlo jet event generator (MAT- 
TER++). Unlike all other previous event genera- 
tors, the current event generator has a space-time 
structure built into the formalism. Similar to Q- 
PYTHIA (as well as JEWEL to some extent) and 
vacuum event generators, but different from MAR- 
TINI, MATTER++ generates events by sampling a 
Sudakov form factor, which is constructed by a com- 
bination of vacuum and medium modified splitting 
functions. Calculations have been carried using both 
a truncated Guo-Wang kernel or the first positive cy- 
cle of the Aurenche-Zakharov-Zaraket kernel. These 
yield similar results, if the transport coefficient is 
modulated with the choice of kernel. Thus the cur- 
rent event generator keeps track of the virtuality of 
every part of the shower explicitly. 

The light-cone locations [x~ locations) of all the 
splits are determined by sampling the Fourier trans- 
form of the uncertainly in the opposite light cone 
momentum (i.e. Ap + ). The width of the x~ dis- 
tribution is set to reproduce the formation time of 
a radiation with a given p~ and p + . As a result, 
the current event generator keeps track of the ex- 
act locations of all the splits in a given event and 
assuming light like propagation between splits, also 
can yield the exact light-cone location of the various 
partons. Transverse broadening is carried out par- 
ton by parton in the regions between the splits, with 
appropriate color factors depending on the flavor of 
the parton. So far calculations have been carried out 
assuming two undistinguished flavors of quark. 

In order to simply study the behavior of the jet 
event generator separated from the dynamics of the 
medium, we have carried out simulations on a static 
medium of finite length, held at a fixed temperature. 
The jets are created at one end and propagate to the 
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other end, showering in the process. Given the lack 
of any transverse dynamics in the medium and the 
desire to not use an jet clustering algorithms, we 
have not discussed the transverse structure of the 
hard jet. The focus has been primarily on the lon- 
gitudinal structure, which may be compared with 
standard DGLAP evolution routines. Using a 6- 
function as the input parton-to-parton fragmenta- 
tion function at the scale of Qo = 1 GeV 2 as well 
as a regular BKK parton-to-hadron fragmentation 
function, we compute and compare the longitudinal 
distributions of quarks and gluons from DGLAP cal- 
culations and the MATTER++ simulations. 

The comparisons for the case of the vacuum in 
Figs.|3]and|4j which represent the case of an originat- 
ing quark jet and a gluon jet respectively, are excel- 
lent. The comparisons with the case of the medium 
with a length of 5 fm and a q = 1 GeV 2 /fm, in 
Figs. [5] and [6j are satisfactory, though not excellent. 
This is due to the somewhat different dynamics that 
is naturally included in a Monte-Carlo simulation 
versus a DGLAP evolution calculation. For exam- 
ple, evolving a (5-function using the DGLAP equa- 
tion is somewhat difficult due to convergence issues. 
Also the full three dimensional integration that is re- 
quired by Eq. (16), can only be carried out approx- 



imately by averaging over the distance dependence 
of the fragmentation function. 

Comparisons using the BKK parton-to-hadron 
fragmentation function show better agreement be- 
tween Monte-Carlo and DGLAP calculations. Even 
this comparison is not exact as partons in the Monte- 
Carlo simulation which drop below the minimum 
Qo — 1 GeV 2 , while more that a fermi away from the 
exit surface are excised from the shower. These can- 
not be fragmented using a perturbatively factorized 
fragmentation function. Also at this low scale these 
interact strongly with the medium and based on 
AdS/CFT calculations, these should be swiftly at- 
tenuated in the medium. We also presented a study 



motivated by the recent measurement of the frag- 
mentation functions within a jet as measured by the 
ATLAS and CMS collaborations. We have clearly 
demonstrated that if some part of the jet's energy is 
lost in the medium and as a result is not present 
in the reconstructed jet energy, then the ratio of 
the medium modified to the vacuum fragmentation 
function shows a typical non-monotonic behavior of 
first decreasing with increasing z and then increasing 
with with increasing z. 

The Monte-Carlo simulations presented in this pa- 
per, are rather simplistic. In future efforts, the 
transverse structure of the jet will be studied us- 
ing jet reconstruction algorithms. The propagation 
of these jets will then be studied in a dynamically 
evolved medium, both in smooth and event-by-event 
fluid dynamical simulations. The production cross 
section of the jets will be calculated from the prod- 
uct of nuclear parton distribution functions and hard 
scattering matrix elements. Finally, a more phe- 
nomenologically accurate hadronization scheme will 
need to be used to hadronize these jets. 
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